function qs = initialize_targeted_subsidies_4(mu,B,A,P,rho,lambda_rnd,lambda_tech)

n = length(A);

%% Computation of the M and R matrices.
M = inv(eye(n)+rho*B-lambda_rnd*A-lambda_tech*P);
q_bar = M*mu;
q_bar = q_bar.*(q_bar>0); % Keep only non-negative values.
R = M*(eye(n)+lambda_rnd*A+lambda_tech*P);

%% Check positive definiteness of the Hessian.
H = eye(n)+2*(eye(n)-1/2*rho*R'*B-R')*R;

%% Compute the analytic targeted subsidies.
subsidy = 2*inv(H+H')*(2*R'*(eye(n)+1/2*rho*B)-eye(n))*q_bar;
subsidy = subsidy.*(subsidy>0);

q_sub = q_bar + R*subsidy;

%% Check if output levels are negative in the counterfactual.
ind = find(q_sub<0);
if(~isempty(ind))
    node_id = [1:n];
    A_prime = A;
    P_prime = P;
    d_prime = sum(A_prime)';
    B_prime = B;
    mu_prime = mu;
    while(~isempty(ind))
        
        %% Remove agents with negative output levels.
        node_id(ind) = [];
        A_prime(ind,:) = [];
        A_prime(:,ind) = [];
        d_prime(ind) = [];
        n_prime = length(A_prime);
        u_prime = ones(n_prime,1);
        P_prime(ind,:) = [];
        P_prime(:,ind) = [];
        B_prime(ind,:) = [];
        B_prime(:,ind) = [];
        mu_prime(ind) = [];
        n_prime = length(A_prime);
        
        %% Computation of the M and R matrices.
        M_prime = inv(eye(n_prime)+rho*B_prime-lambda_rnd*A_prime-lambda_tech*P_prime);
        q_bar_prime = M_prime*mu_prime;
        q_bar_prime = q_bar_prime.*(q_bar_prime>0); % Keep only non-negative values.
        R_prime = M_prime*(eye(n_prime)+lambda_rnd*A_prime+lambda_tech*P_prime);
        
        %% Check positive definiteness of the Hessian.
        H_prime = eye(n_prime)+2*(eye(n_prime)-1/2*rho*R_prime'*B_prime-R_prime')*R_prime;
        
        %% Compute the analytic targeted subsidies.
        subsidy_prime = 2*inv(H_prime+H_prime')*(2*R_prime'*(eye(n_prime)+1/2*rho*B_prime)-eye(n_prime))*q_bar_prime;
        subsidy_prime = subsidy_prime.*(subsidy_prime>0);
        
        q_sub_prime = q_bar_prime + R_prime*subsidy_prime;

        ind = find(q_sub_prime<0);
        
    end
    subsidy = zeros(n,1);
    q_sub = zeros(n,1);
    for i=1:length(node_id)
        subsidy(node_id(i)) = subsidy_prime(i);
        q_sub(node_id(i)) = q_sub_prime(i);
    end
    
end

qs = [q_sub; subsidy];
